clearvars -except do_QC

date     = SAT_SST_func_IO('date');
if do_QC == 0
    P.app      = 'qcu';
else
    P.app      = 'qcf';
end

% *************************************************************************
% Read and parse temperature 
% *************************************************************************
dir      = SAT_SST_func_IO('GHCN');
dir_save = SAT_SST_func_IO('GHCN_processed',P);

file = [dir,'ghcnm.tavg.v4.0.1.',date,'.',P.app,'.dat'];

fid = fopen(file,'r');
ch = fscanf(fid,'%c',10000000000000);
fclose(fid);

data = reshape(ch,116,numel(ch)/116)';
data = data(:,1:end-1);
station = data(:,1:11);
yr   = str2num(data(:,12:15));
var  = data(:,16:19);

value = nan(size(data,1),12);
for mon = 1:12
    temp = data(:,[20:24]+(mon-1)*8);
    value(:,mon) = str2num(temp)/100;
end
value(value == -99.99) = nan;

[uni_station,~,J] = unique(station,'rows');

% *************************************************************************
% Read and parse metadata 
% *************************************************************************
file_inv = [dir,'ghcnm.tavg.v4.0.1.',date,'.',P.app,'.inv'];
fid = fopen(file_inv,'r');
ch = fscanf(fid,'%c',10000000000000);
fclose(fid);

clear('inv')
list = [0 find(ch == 10)];
for ct = 1:numel(list)-1
    temp = ch(list(ct)+1:list(ct+1));
    inv(ct,:) = temp(1:68);
end
inv_ID = inv(:,1:11);
lat = str2num(inv(:,13:20));
lon = str2num(inv(:,22:30));
lev = str2num(inv(:,32:37));
lev(lev > 9900) = nan;

% *************************************************************************
% save data
% *************************************************************************
for ct = 1:size(uni_station,1)
    
    try
        clear('data_yr','data_T','l','data_lon','data_lat','data_lev')
        l = J == ct;
        data_yr = yr(l);
        data_T  = value(l,:)';


        l = find(ismember(inv_ID,uni_station(ct,:),'rows'));
        data_lon = lon(l);
        data_lat = lat(l);
        data_lev = lev(l);

        file_save = [dir_save,uni_station(ct,:),'_mon_V4.mat'];
        save(file_save,'data_yr','data_T','data_lon','data_lat','data_lev','-v7.3')
    catch
        disp(['something is wrong for station ', uni_station(ct,:)])
    end
end
